histDist1<-function (y, family = NO, freq = NULL, density = FALSE, nbins = 10, 
    xlim = NULL, ylim = NULL, main = NULL, xlab = NULL, ylab = NULL, 
    data = NULL, ...) 
{
    FA <- as.gamlss.family(family)
    fname <- FA$family[1]
    dfun <- paste("d", fname, sep = "")
    lpar <- length(FA$parameters)
    typeDist <- FA$type
    subsY <- substitute(y)
    y <- if (!is.null(data)) 
        get(deparse(substitute(y)), envir = as.environment(data))
    else y
    freq <- if (!is.null(data) && !is.null(freq)) 
        get(deparse(substitute(freq)), envir = as.environment(data))
    else freq
    switch(typeDist, Continuous = {
        extra <- (max(y, na.rm = TRUE) - min(y, na.rm = TRUE))/5
        xmin <- if (is.null(xlim)) min(y, na.rm = TRUE) - extra else xlim[1]
        xmax <- if (is.null(xlim)) max(y, na.rm = TRUE) + extra else xlim[2]
        if (!FA$y.valid(xmin)) xmin <- 0.001
        if (!FA$y.valid(xmax)) xmax <- 0.999
        x1 <- seq(xmin, xmax, length = 101)
    }, Discrete = {
        if (fname %in% .gamlss.bi.list) {
            if (NCOL(y) == 1) {
                x1 <- c(0, 1)
                bd <- rep(1, 2)
            } else {
                if (any(abs(y - round(y)) > 0.001)) {
                  warning("non-integer counts in a binomial GAMLSS!")
                }
                bd <- y[, 1] + y[, 2]
                y1 <- y[, 1]
                if (any(y1 < 0 | y1 > bd)) stop("y values must be 0 <= y <= N")
                xmin <- if (is.null(xlim[1])) {
                  if (all(bd == bd[1])) min(y1, na.rm = TRUE) else 0
                } else xlim[1]
                xmax <- if (is.null(xlim)) {
                  if (all(bd == bd[1])) max(bd, na.rm = TRUE) else 1
                } else xlim[2]
                x1 <- seq(xmin, xmax, by = 1)
            }
        } else {
            xmin <- if (is.null(xlim)) min(y, na.rm = TRUE) else xlim[1]
            xmax <- if (is.null(xlim)) max(y, na.rm = TRUE) else xlim[2]
            x1 <- seq(xmin, xmax, by = 1)
        }
    }, Mixed = {
        stop("Mixed distributions are not implemented yet")
    })
    if (is.null(freq)) {
        {
            mod <- try(gamlssML(y, family = fname), silent = TRUE)
            if (any(class(mod) %in% "try-error")) {
                mod <- try(gamlss(y ~ 1, family = fname, ...))
            }
        }
    }
    else {
        {
            mod <- try(gamlssML(y, weights = freq, family = fname, 
                ...), silent = TRUE)
            if (any(class(mod) %in% "try-error")) {
                mod <- try(gamlss(y ~ 1, weights = freq, family = fname, 
                  ...))
            }
        }
    }
    mod$call$family <- eval(as.expression(fname))
    if (mod$method == "BFGS" || mod$method == "nlminb") {
        mod$call$formula <- subsY
    }
    else {
        mod$call$formula[[2]] <- subsY
    }
    if (!is.null(data)) 
        mod$call$data <- substitute(data)
    switch(lpar, `1` = {
        newcall <- if ((fname %in% .gamlss.bi.list)) call(dfun, 
            x1, mu = fitted(mod)[1], bd = bd[1]) else call(dfun, 
            x1, mu = fitted(mod)[1])
    }, `2` = {
        newcall <- if ((fname %in% .gamlss.bi.list)) {
            call(dfun, x1, mu = fitted(mod)[1], sigma = fitted(mod, 
                "sigma")[1], bd = bd[1])
        } else {
            call(dfun, x1, mu = fitted(mod)[1], sigma = fitted(mod, 
                "sigma")[1])
        }
    }, `3` = {
        newcall <- if ((fname %in% .gamlss.bi.list)) {
            call(dfun, x1, mu = fitted(mod)[1], sigma = fitted(mod, 
                "sigma")[1], nu = fitted(mod, "nu")[1], bd = bd[1])
        } else {
            call(dfun, x1, mu = fitted(mod)[1], sigma = fitted(mod, 
                "sigma")[1], nu = fitted(mod, "nu")[1])
        }
    }, `4` = {
        newcall <- if ((fname %in% .gamlss.bi.list)) {
            call(dfun, x1, mu = fitted(mod)[1], sigma = fitted(mod, 
                "sigma")[1], nu = fitted(mod, "nu")[1], tau = fitted(mod, 
                "tau")[1], bd = bd[1])
        } else {
            call(dfun, x1, mu = fitted(mod)[1], sigma = fitted(mod, 
                "sigma")[1], nu = fitted(mod, "nu")[1], tau = fitted(mod, 
                "tau")[1])
        }
    })
    switch(typeDist, Continuous = {
        y1 <- eval(newcall)
        xlim <- c(xmin, xmax)
        main <- if (is.null(main)) paste("The ", deparse(subsY), 
            " and the fitted ", FA$family[1], " distribution", 
            sep = "") else main
        if (!is.null(freq)) y <- rep(y, freq)
        xlab <- if (is.null(xlab)) deparse(subsY) else xlab
        ylab <- if (is.null(ylab)) paste("f()") else ylab
        truehist(y, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, 
            nbins = nbins, main = main, col = "gray", lty = 3, 
            border = "blue", fg = rainbow(12)[9], col.main = "blue4", 
            col.lab = "blue4", col.axis = "blue")
        lines(x1, y1, col = "red", col.axis = "blue", col.main = "blue4", 
            col.lab = "blue4", lwd = 2, fg = gray(0.7))
        if (density == TRUE) {
            dens <- density(y)
            lines(dens$x, dens$y, col = "blue")
        }
    }, Discrete = {
        if (fname %in% .gamlss.bi.list) {
            if (all(bd == bd[1])) {
                y1 <- eval(newcall)
                tabley <- if (NCOL(y) == 1) table(y) else table(y[, 
                  1])
                dft <- data.frame(tabley)
                r <- barplot(tabley/sum(tabley), fg = "blue", 
                  col = "gray", axis.lty = 1, border = "blue", 
                  col.axis = "blue4", col.main = "blue4", col.lab = "blue4", 
                  ylim = ylim, xlim = xlim, ylab = ylab, xlab = xlab, 
                  main = main)
                yy1 <- y1[x1 %in% dft[[1]]]
                lines(r, yy1, type = "h", col = "red", lwd = 2, 
                  col.axis = "blue", col.main = "blue4", col.lab = "blue4", 
                  fg = gray(0.7))
                points(r, yy1, col = "red", pch = 21, lwd = 2, 
                  col.axis = "blue")
            } else {
                xlim <- c(xmin, xmax)
                main <- paste("proportions of ", deparse(subsY), 
                  sepp = "")
                main <- if (is.null(main)) paste("proportions of ", 
                  deparse(subsY), sepp = "") else main
                xlab <- if (is.null(xlab)) deparse(subsY) else xlab
                ylab <- if (is.null(ylab)) paste("f()") else ylab
                truehist(y1/bd, xlim = xlim, xlab = "proportions", 
                  ylab = ylab, ylim = ylim, nbins = nbins, col = "gray", 
                  lty = 3, border = "blue", col.axis = "blue", 
                  col.main = "blue4", col.lab = "blue4", main = main, 
                  fg = rainbow(12)[9])
                lines(fitted(mod), rep(1, length(fitted(mod))), 
                  type = "h", col = "red")
                points(fitted(mod), rep(0, length(fitted(mod))), 
                  col = "red", pch = 25)
                points(fitted(mod), rep(1, length(fitted(mod))), 
                  col = "red", pch = 24)
            }
        } else {
            y1 <- eval(newcall)
            fy <- if (is.null(freq)) factor(y, levels = x1) else factor(rep(y, 
                freq), levels = x1)
            notresy <- if (is.null(freq)) factor(y) else factor(rep(y, 
                freq))
            dft <- data.frame(tabley <- xtabs(~fy))
            main <- if (is.null(main)) paste("Barplot of the ", 
                deparse(subsY), " and the fitted ", FA$family[2], 
                " distribution", sep = "") else main
            r <- barplot(tabley/sum(xtabs(~notresy)), fg = "blue", 
                col = "gray", axis.lty = 1, border = "blue", 
                col.axis = "blue4", col.main = "blue4", col.lab = "blue4", 
                main = main, ylim = ylim, ylab = ylab, xlab = xlab)
            yy1 <- y1[x1 %in% dft[[1]]]
            lines(r, yy1, type = "h", col = "red", lwd = 2, col.axis = "blue", 
                col.main = "blue4", col.lab = "blue4", fg = gray(0.7))
            #points(r, yy1, col = "red", pch = 21, lwd = 2, col.axis = "blue")
        }
    }, Mixed = {
        stop("Mixed distributions are not implemented yet")
    })
    mod$yy1<-yy1
mod$tt1<-tabley/sum(xtabs(~notresy))
mod
}